Loop updates for variational and projector quantum Monte Carlo simulations 

in the valence-bond basis 
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We show how efficient loop updates, originally developed for Monte Carlo simulations of quantum 
spin systems at finite temperature, can be combined with a ground-state projector scheme and vari- 
ational calculations in the valence bond basis. The methods are formulated in a combined space of 
spin z-components and valence bonds. Compared to schemes formulated purely in the valence bond 
basis, the computational effort is reduced from up to 0(N 2 ) to O(N) for variational calculations, 
where N is the system size, and from 0(m 2 ) to 0(m) for projector simulations, where m 3> N 
is the projection power. These improvements enable access to ground states of significantly larger 
lattices than previously. We demonstrate the efficiency of the approach by calculating the sublattice 
magnetization M s of the two-dimensional Heisenberg model to high precision, using systems with 
up to 256 x 256 spins. Extrapolating the results to the thermodynamic limit gives M s — 0.30743(1). 
We also discuss optimized variational amplitude-product states, which were used as trial states in 
the projector simulations, and compare results of projecting different types of trial states. 

PACS numbers: 02.70. Ss, 75.10.Jm, 75.40.Mg, 75.40.Cx 



I. INTRODUCTION 

An ongoing challenge in simulations of quantum spin 
systems is to reach larger lattices sizes, thus enabling 
more reliable extrapolations to the thermodynamic limit. 
With the advent of loop-cluster algorithms 1-5 and related 
schemes^ developed since the mid-1990s, finite temper- 
ature (T) quantum Monte Carlo (QMC) simulations have 
become possible on lattices with millions of spins for 
models with positive-definite path integral (world line)&2 
or stochastic series expansion (SSE) 1 - representation of 
the partition function. The computational effort scales 
linearly in the number of spins N. Since the effort also 
scales as 1/T, simulations at very low T or in the ground 
state (using T low enough to eliminate finite-T effects), 
are limited to smaller lattices, however. The ground state 
can typically be reached for « 10 4 spins. 

Currently accessible system sizes suffice for studying 
ground states of many important models, e.g., the two- 
dimensional (2D) Heisenberg mode l 1 1 i 12 and variants of 
it with non-uniform coupling patterns leading to quan- 
tum phase transitions of the antiferromagnet into a dis- 
ordered ground state i 13 ' 14 In other, similar systems there 
are still controversial issue a 15 ' 16 that may need larger lat- 
tices to be conclusively resolved. Larger lattice sizes are 
crucial in systems exhibiting more complex ground states 
and quantum phase transitions. One example of cur- 
rent interest is a class of "J-Q" models — 2D Heisenberg 
models with four-spin interactions engineered to destroy 
the antiferromagnetic order and drive the system into a 
valence-bond-solid (VBS) statei£"i^ The VBS state has 
intricate fluctuations and the true nature of its thermo- 
dynamic limit is only manifested on large lattices^ This 
illustrates the need to develop better ground-state meth- 
ods, as a more efficient alternative to going to very low 
T with finite-temperature methods. 



In this paper we introduce a method combining loop 
updates first developed for finite-T simulations 1 ^' 5 - with 
a ground-state projector QMC method operating in the 
valence bond (VB) basis i 20 ' 21 This over-complete sin- 
glet basis has some features that make it uniquely well 
suited for studies of spin-rotationally invariant hamilto- 
nians such as the Heisenberg model and its extensions 
with multi-spin interactions *2. 

It has been known for some time^ that there is a simple 
and elegant relationship between VB states consisting of 
N/2 pairs of spins forming singlets^ 2 , and the loop algo- 
rithm, which indeed works by switching between a VB 
basis and a basis of N spins f and J, (for S = h systems). 
Here we exploit this switching for ground state projec- 
tions. An attractive feature of this approach is that it 
enables the use of very good singlet trial wave functions, 
the simplest example of which is the amplitude-product 
state proposed by Liang, Doucot and Anderson i 23 ' 24 The 
ground state can then be reached much faster than with 
finite-T methods, and with much less computational ef- 
fort than projector methods formulated purely in the VB 
basisi 20 ' 25 ' 26 In addition we show that purely variational 
calculations can also be made more efficient by combin- 
ing spins and VBs, including a loop update similar to one 
previously developed for classical dimer modelsj 28 i 29 

The projector QMC algorithm with loops, working in 
the combined space of VBs and spins, is in the end very 
similar to T > SSE and worldline loop algorithms. 
Essentially, the T = projector approach corresponds 
to "cutting open" the periodic imaginary-time boundary 
and "sealing" the resulting open loop segments with va- 
lence bonds (which serve as continuations of the loops). 

We demonstrate the efficiency of the projector method 
by producing high-precision bench-mark results for the 
sublattice magnetization of the 2D Heisenberg model 
with up to 256 x 256 spins. We also discuss the prop- 



erties of the variationally optimized amplitude-product 
states used as a trial states for the ground-state projec- 
tions (extending the results of Ref. [24] to larger lattices) . 
We begin in Sec. [TT] by summarizing the properties 
of the VB basis needed for formulating the algorithms. 
In Sec. IIIII we discuss variational Monte Carlo simula- 
tions and optimization of amplitude-product states, and 
in Sec. II VI we describe the projector QMC method. We 
show results of both variational and projector calcula- 
tions in Sec. El and conclude in Sec. |VI] with a summary 
and discussion. 



II. THE VALENCE BOND BASIS 



A VB state is a product of two-spin singlets, 



M) = (| t«-U> " Uati>)/V2, 



(1) 



where a and b denote sites on sublattice A and B on a 
bipartite system. For N (even) spins there are (AT/2)! 
ways to draw the VBs, and hence the basis is massively 
overcomplete. The properties of VB states have been 
discussed extensively in the literature ! 23 ' 32 ' 33 Here we will 
only summarize the most important aspects of the basis 
and introduce some notation needed for our algorithms. 
We can formally use a label r = 1, . . . , (N/2)\ for enu- 
merating the VB configurations and denote a state as 



\V r 



\{al,bl){al,bl)---{ar m ,V m )). 



The overlap between two VB states is 

{V l \V r ) = 2 N ^- N '\ 



(2) 



(3) 



where JVi oop is the number of loops formed in the so-called 
transposition graph when the bonds in |Vj) and \V r ) are 
superimposed. An example is shown in Fig. [T] 

Like the overlap, matrix elements of operators of in- 
terest can typically also be related to the transposition- 
graph loops i 23 ' 32 ' 33 To compute spin-spin correlations we 
will need 



(Vl\V r ) 



0. 



if Xi ^ Xj , 



0^(3/4), if A, -A, 



(4) 



where <fn = ±1 for sites on sublattice A and B, respec- 
tively, and Aj € {l,Ai oop } is a label distinguishing the 
loop to which site i belongs. More complicated matrix 
elements and their relationships to the loop structure are 
discussed in Ref. |33|. 

Fig. Q] also shows one of the spin states, 

\ZT) = \S!(r,i),...,S* N (r,i)), i = 1,.. . ,2 N '*, (5) 

contributing to both the VB states (meaning that the 
spins on all bonds of the transposition graph are anti- 
parallel). In ((SJ) the index i refers to the allowed spin 
states in \V r ); the 2 N / 2 states with anti-parallel spins on 




FIG. 1: (Color online) Transposition graph of VB states on a 
5x4 lattice with sublattices A and B indicated with open and 
solid circles. Bond configurations of a bra (Vi\ and ket state 
\V r ) are shown as dashed and solid lines, respectively. For 
clarity, all bonds here are of length one lattice spacing, but 
in general a bond can connect any pair of sites on different 
sublattices. The number of loops is N\ oop = 3 and the number 
of sites N = 20, giving, according to Eq. (3J, an overlap 
(FilVv) = 2~ 7 . A configuration of f and 4- spins compatible 
with both VB states is also shown. 

every bond. With the sign convention in (JXJ) , the VB 

state can be written as 



\Vr) 



1 



2 N/4 



2 JV/2 

E(- 



■ l) A t( r A\Zl 



(6) 



where Af(r, i) denotes the number off spins on sublattice 
A. The VB states thus obey the well known Marshall sign 
rule, which ensures that the ground state has a positive- 
definite expansion in terms of them [for an interaction 
with the same A-B bipartite structure as in the definition 
of the valence bonds, Eq. (JTJ]. 

When writing the overlap Q using spin states, 



{Vi\V r ) = 



1 



2 JV/2 



^(ZJ|Z[)(-l) A t(^)+^('j) ; 



(7) 



S3 



only the terms with Z\ = Z= contribute. Since the spins 
on every bond must be anti-parallel, the spin configura- 
tion around a loop in the overlap graph must be stag- 
gered. The signs cancel and the sum then simply counts 
the number of such configurations. There are two stag- 
gered spin patterns for each loop, resulting in the overlap 
([3]). Here we will use the equivalence between the two 
ways of expressing the overlap, © and J7]), as a starting 
point for constructing efficient variational and ground- 
state projector algorithms. 



III. VARIATIONAL MONTE CARLO 

We first discuss variational calculations, which we will 
use to construct good trial states for subsequent ground- 
state projection. An arbitrary singlet state |\E') can be 
expanded in VB states; 



i*>=5~> r |v P ). 



(8) 



Because of the over-completeness, the expansion coeffi- 
cients w r are not unique, but this is not a problem in 
practice. In the variational state introduced by Liang et 
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FIG. 2: (Color online) A two-bond update. For any pair of 
two sites on the same sublattice (here the sites a, c, as in- 
dicated by the brackets), the two bonds connected to them 
can be reconfigured in a unique way while maintaining the re- 
striction of bipartite bonds. Detailed balance is then satisfied 
for updates toggling between the two configurations, with the 
Metropolis acceptance probability (|12p . 

al.jQ the coefficients are taken to be products of bond 
amplitudes h(r) > 0, 



n^r) 



N r (r 



(9) 



where N r (r) is the number of bonds of size r, where in 
a translationally-invariant system, all applicable lattice 
symmetries can be used, e.g., for a 2D square lattice the 
amplitudes are h(x,y), with x = \r x \ and y — \r y \ (the 
components of the vector r defining the "shape" of the 
bond, transformed to the all-positive quadrant), and also 
h(y,x) = h{x,y). 

To optimize the amplitudes using variational QMC 
simulations, the energy is written as 



E 



ElrWlr 



(tfltf) 



(10) 



where the weight Wi r and energy estimator E\ r are 

(Vi\H\V r ) 



Wl r = WlW r {Vl\V r ) , E lr = 



(Vl\V r ) 



(11) 



For a model with two-spin interactions (to which the 
methods discussed here are not limited ) 17 ' 20 but for sim- 
plicity of the discussion we will not consider higher-order 
interactions here) , the energy estimator is a sum of terms 
of the form (|1J . As in any variational calculation, the idea 
is to compute E and its derivatives with respect to the 
amplitudes and then use this information to periodically 
adjust the amplitudes, in order to approach the energy 
minimum. Here we will not be concerned with the op- 
timization aspect of the problem (which is discussed in 
detail in Ref. [24], including also explicit expressions for 
evaluating the derivatives) , but focus on the Monte Carlo 
sampling aspects of the problem. 



A. Monte Carlo sampling 

For a given set of amplitudes, the VB configurations 
can be sampled according to the weight Wj r by simple 
moves of bond pairs £2- Choosing two sites on the same 
sublattice (e.g., next- nearest- neighbors on the square lat- 
tice), there are two allowed (bipartite) configurations of 
the two bonds connected to these sites, and the Monte 
Carlo move amounts to switching from the current one to 



the other possible one. Such an update is illustrated in 
Fig. H The ratio W Vr >/Wi r of the weights after (W/v) 
and before (W; r ) the update needed for the Metropolis 
acceptance probability; 



P 



accept 



W v 



Ww 



,1 



(12) 



The ratio involves the amplitudes of the affected bonds 
from the amplitude product ((5J, as well as the ratio of 
the new to old overlap ((3J). Here the number of loops 
iVioop can increase or decrease by one. Determining this 
change requires tracing a loop going through one of the 
four sites involved in the two-bond move. Starting from 
one of these sites, the loop going through that site before 
the re-configuration is followed, until some other site out 
of those four is reached. If that site is connected to the 
same bond as the initial site, then the four sites must 
belong to two different loops (and after the move they 
will all be in the same loop, i.e., Afi 00 p is decreased by 
1), whereas if it belongs to the second bond all four sites 
are in a single loop (which is split into two after the up- 
date, leading to N\ oop increasing by 1). The loop tracing 
is the most time consuming part of the calculation, in 
particular for systems where the loops are typically long. 
The loop lengths are related to the squared sublattice 
magnetization, according to 32 
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(13) 



A=l 



where L\ is the length of loop A. For the single loop of 
length Li going through a randomly chosen site i, this 
relation becomes 2 -^ 



(L i ) = N(M%) 



(14) 



where the averages refer to the Monte Carlo sampling. 
For systems with antiferromagnetic long-range order 
there are system-spanning loops; 23 ' 32 i.e., there are some 
loops of length L\ ex A in typical configurations. In this 
case, the computational effort of a sweep of oc A" two- 
bond updates therefore scales as A' 2 . 

Here we avoid the loop-tracing step by re-introducing 
the spins into the sampling space, replacing the pure- 
bond overlap © with the equivalent expression (|7J|. In 
addition to sampling the bond tilings Vi , V r of the bra and 
ket states, we then also have to sample their compatible 
spin configurations Z\ = Z- = Z* r . A two-bond move 
cannot always be performed in this combined space, since 
the re-configured bonds may not be compatible with the 
current spin configuration. Starting with an allowed com- 
bined configuration (Vi, V r , Z\ r ) (an example of which is 
shown in Fig. [1]), we carry out N/2 two-bond attempts 
in each of the bond configurations Vi,V r , and thereafter 
construct all VB loops and flip each loop (i.e., flip all 
the spins in the loop) with probability 1/2. The number 
of operations required for one such full cycle of updates 
scales as N. In the scaling of the computational effort we 



thus avoid the factor (L\ oop ) = N(M^). If there is anti- 
ferromagnetic order (system spanning- loops, L\ oop ex N), 
or even in the absence of long-range order if the correla- 
tion length is large, the time savings of switching to the 
combined spin-bond basis is very large. 

It is worth pointing out the reason why this method of 
avoiding to compute the overlap ([3]) works. It is because 
the combined Monte Carlo sampling of spins and bonds 
will, according to the standard detailed balance and er- 
godicity principles, automatically favor bond configura- 
tions for which there are more compatible spin configura- 
tions, on average exactly in proportion to the re- written 
overlap, Eq. 0, which is just the number of spin con- 
figurations compatible with a given bond configuration. 
Thus, the overlap is taken into account probabilistically 
in the extended space, instead of being computed exactly 
in the pure VB spaced One may ask whether this could 
lead to a worse performance of the combined bond-spin 
sampling relative to pure VB sampling. The bond sam- 
pling in the mixed space is constrained, but, on the other 
hand, in the pure VB sampling a decreasing overlap af- 
ter a bond reconfiguration (which should occur in half 
of the updates) reduces the acceptance rate. These two 
effects should, on average, balance each other, and thus 
the performance of the two methods should be similar 
when measured in terms of the number of bond updates 
performed. The sampling is, however, much faster in the 
combined space. 

Note also that measurements of observables can (and 
normally should) still be carried out in the pure VB ba- 
sis (instead of measuring just the z-components of, e.g., 
correlation functions), i.e., at this stage it does not mat- 
ter what the spin configuration is, and one just considers 
rotationally invariant estimators, such as Eq. ((4]), in the 
pure VB basis. This corresponds exactly to summing 
over all compatible spin configurations, i.e., one can con- 
sider the VB (loop) estimators as improved estimators^ 



B. Bond-loop updates 

As an alternative to sampling the VB states using 
the two-bond reconfigurations, one can also implement a 
loop update similar to one developed for classical dimer 
modelSi 28 i 29 The idea is to first move one end of a ran- 
domly chosen bond, thereby creating two defects (an 
empty site and one site with two bonds) that should not 
be present in a valence-bond state. Subsequent bond 
moves are then carried out to move one of these defects, 
until it annihilates the second one (the loop closes) and a 
new allowed bond configuration has been generated. The 
first step of such a loop update is illustrated in Fig. [3] 
(Note that these loops are different from those in the rest 
of the paper). 

In an algorithm (as well as in an actual computer im- 
plementation), it is convenient to represent the bonds 
by links between sites, stored in a one-dimensional array 
v(i), i = 1,...,N. Thus, if sites i and j are connected 
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FIG. 3: (Color online) The first step in a bond-loop update. 
Open and solid circles indicate up and down spins compati- 
ble with the valence bonds. Given an initial site jo, the site 
i linked to it is identified. The bond jo,i is then replaced 
by i,j, with j on the opposite sublattice from i and chosen 
probabilistically as discussed in the text. If Si = Sj, the 
new bond configuration is not allowed, and a new j should be 
generated. The bond move results in one site with no bond 
and one with two bonds (unless j — jo and the loop update 
is completed). The site linked by the old bond to the site 
with two bonds becomes site i for the following step. These 
procedures are repeated until the loop closes (j —jo). 



by a bond, then v(i) = j and v(j) = i. Starting at 
a randomly chosen site jo (in either \Vi) or \V r ); here 
we consider \V r ) with the links stored in an array v r ) 
the bond connected to it is considered for a reconfigu- 
ration. The bond will stay attached to its second site, 
i — v r (jo), and the new other end, j, will be chosen ac- 
cording to probabilities proportional to the correspond- 
ing amplitudes h{vij). This is accomplished in practice 
by carrying out a bisection search for values bracketing a 
random number in [0,1) in a pregenerated table contain- 
ing cumulative probabilities, i.e., normalized partial sums 
of the amplitudes^ If, for the chosen site j [obtained by 
adding the chosen bond vector r to the current site i] the 
spin S z (j) = S z {i), the move would not be consistent 
with the spin configuration, and a new j should then be 
chosen [repeatedly if necessary, until S z (j) ^ S z (i) is sat- 
isfied]. After moving the bond to an acceptable j, setting 
v r {i) = j, this site will have two bonds on it (but in the 
array v r the original link has been destroyed, which does 
not matter since it will no longer be needed), and the 
original site jo will have no bond attached to it (unless 
the chosen site happens to be the starting site, j = jo, 
in which case the loop update is immediately terminated 
and another loop building is started from a new randomly 
chosen site 3 -). To "heal" the double-bond defect at j, the 
end of the old bond connected to it is moved. This is done 
in the same way as for the initial bond move, with jo in 
the description above replaced by j. The procedures are 
repeated until it happens that j = jo , at which point the 
double-bond and no-bond defects annihilate each other 
and the loop closes. The loops can be large, and this 
kind of update should therefore be much more efficient 
than the local two-bond updates. The number of loops 
constructed in each updating sweep should be adjusted 
so that, on average, a number oc N of bonds are moved. 

For an amplitude-product state with non-zero ampli- 



tildes /i(r) for all r, the bisection search in the bond-loop 
update introduces a factor oc ln(7V) in the time scaling of 
the algorithm, which may or may not (depending on the 
nature of the state) be outweighed by shortened auto- 
correlation times relative to the two-bond updates. Note 
also that if the amplitudes are zero or very small beyond 
some distance r (e.g., in the case of exponentially decay- 
ing amplitudes), the table of probabilities should only 
include those amplitudes with practically non-vanishing 
amplitudes, which also removes the ln(iV) factor from the 
scaling in such cases. 



IV. PROJECTOR MONTE CARLO 

We now turn to the projector QMC method, using the 
Heisenberg model as a concrete example of loop updates 
for efficient ground-state projection of a variational state 
in the VB basis. We write the hamiltonian on a bipartite 
lattice as 



H — — 2__, H a b 

(a,b) 



H a b — — (S a • Si, — -7), 



(15) 



where (a, b) denotes nearest-neighbor sites. The bond 
operators H a b are singlet projectors (equivalent to loop 
operators, 2 up to a factor of 2 and bipartite rotation) 
which can have two effects when acting on a VB state: 

Ha6(o,6) = (o,6), (16) 

ff od (a,&)M) = |(M)(c,6). (17) 

These simple rules, and the absence of minus signs (in the 
case of an unfrustrated system), enable a QMC scheme 
for projecting out the ground state |\&o) from an arbitrary 
state dSJ in the VB basis.^^^ 



A. Ground state projection 

Irrespective of the basis and hamiltonian, the projector 
approach is based on the expression 



(-H) m \Vr) = c (-E y 



|0> 



A-l 



71=1 



C() 



E n 

E 



(18) 

where |n), n — 0, . . . , A — 1 are the energy eigenstates 
in a Hilbert space of A states. If Eq is the largest (in 
magnitude) eigenvalue and the expansion coefficient Cq 7^ 
0, then |0) oc {-H) m \^) for large m. In a QMC projector 
scheme, a high power of H and its action on the trial state 
|\t) are sampled stochastically. 

To this end, for the Heisenberg model we write {—H) m 
as a sum of all products of m bond operators and intro- 
duce the notation P a for such operator strings; 



(-#)"' = £ II #«?*? = E P - 



(19) 



where a is a formal label for the different strings of singlet 
projectors. We write the state resulting when acting with 
a P a on a given VB state \V r ) as 



P a \Vr) = {\)°-\V r {a)), 



(20) 



where |V r (a)) is obtained in practice by successively ap- 
plying the rules ([TBI) and (1X7)) , which also gives the num- 
ber o r a of off-diagonal operations (fT7|) . The expectation 
value of an arbitrary operator O can be written as 



(O) 



(*\{-HrO(-H) m \V) = £ ira/9 W£ 
(y\(-H)*™\*) £ 






where the weight W^ r p and estimator O'j'f are 



Irafi "It 
a/3 



a/3 



Ir 



wiw r (Vi(/3)\V r (a))2- 



(°"+°r 



Of = <VK/3)|0|W)/W)|K(a)). 



(21) 



(22) 
(23) 



The expectation value can be evaluated by importance- 
sampling, as discussed in Ref. l2Ct However, up until 
now the sampling was rather inefficient. Each random 
replacement of operators (one or several; the number is 
adjusted to give a suitable acceptance rate) in P a re- 
quires a full propagation of the current VB state, count- 
ing o r a in Eq. (f20|). and thereafter counting Vi oop in the 
overlap ([3]). That is, unlike what is normally the case 
in Monte Carlo simulations, one here cannot simply de- 
termine the weight ratio locally when a small change is 
made in the configuration, and, hence, the full weight has 
to be computed in each update. This results in a scaling 
~ max(m 2 , Nm) of the computational effort involved in 
a full updating sweep. Here one factor m corresponds to 
the effort of propagating the full updated operator string. 
Constructing and counting the loops in the transposition 
graph ([3]) requires oc N operations. If m 3> N this ef- 
fort is swamped by the bond propagation, and the total 
effort of carrying out a number oc m of updates (needed 
to significantly change the operator sequence) is oc m 2 , 
whereas for the situation N > m (which is less interest- 
ing in practice, as we will see when discussing the con- 
vergence properties), the effort is formally oc Nm. When 
using an amplitude-product as the trial state (instead 
of a fixed bond configuration), we should also perform 
oc N two-bond updates of the trial state, each of which 
demands the same computational effort as an operator 
replacement. 



B. Combined bond-spin space 

We will now show how using the combined spin-bond 
basis brings the scaling of one updating sweep down from 
max(m 2 , Nm) to ~ max(V, m). We construct a loop up- 
date, which we here implement for the power (—H) 2m in 
a way similar to the "operator-loop" update for e~@ H , 
where /3 = 1/T, in the SSE representation 5 in T > 
simulations. The operator-loop approach is analogous 




FIG. 4: (Color online) A VB-spin-operator configuration con- 
tributing to (*|(-#) 2m |*) for a 4-site system with m — 2. 
The arcs to the left and right indicate VB states {Vi\, \V r ) and 
the two columns of filled and open circles represent •f and 4- 
spins of compatible spin states (Zj\, \Zj). The spins at the 
four operators (vertices) are also indicated. There are three 
loops, part of which consist of VBs. Expectation values are 
evaluated at the mid-point indicated by the dashed line. 



to the original loop method for world lines^ with the 
main difference that the sampling scheme is formulated 
in terms of insertions and removals of operators, instead 
of deforming world-line configurations. The formal re- 
lationships between the two approaches are discussed in 
Ref. 0- The SSE formulation is somewhat easier to relate 
directly to the present case of projector QMC. We could 
also project with e _/3ff , Taylor-expanded as in the SSE 
method, but the fixed-power approach allows for some 
minor simplifications. The primary differences with re- 
spect to finite- T SSE simulations is then the fixed value 
of m and the VB "end-cap" boundaries in the "propa- 
gation" direction, replacing the periodic boundary con- 
ditions appropriate at finite T. 

To make the analogy with the SSE loop method as 
close as possible, we split the singlet-projector operators 
H a b into their diagonal, H a i(l), and off-diagonal, H a b{2), 
parts, 



#a&(l) = (4 _ s l s b) 
H ab (2) = -i(S'+S , b " 



&a ^b ) 



(24) 
(25) 



We use a superscript e on the operator string P^ in (Til)]) 
to denote the 2'™ different combinations of diagonal and 
off-diagonal operators for a given full-operator string P a . 
With spins Z\, Zj compatible with V r , Vi, we sample 
VBs, spins, and operators, according to the weight 



TV 



«/3,e/ _ 



mw r a) 2m+N/2 oc Wl t 



(26) 



under the condition P^\ZT) — Pl\ZA. The constraints 
exactly compensate for the other factors in the original 
weight ([22|) and there is no explicit dependence in (f26| on 
the operator-string (a,/3, e, /) and spin (i,j) indices. An 
example configuration is shown in Fig. 2J On a bipartite 
lattice the weights are positive, since minus signs present 
in the states © compensate those arising from an odd 
number of off-diagonal operators (1251) (or, equivalently, 
all signs could be eliminated by a sublattice rotation 2 ). 



C. Sampling method 

We now briefly describe the Monte Carlo sampling pro- 
cedures. Starting with VB configurations V r ,Vi (where 
normally one would take V r — VJ for simplicity) and com- 
patible spin configurations Z r = Z , an initial string co- 
taining only diagonal operators H a b(l) can be used (con- 
sistent with the constraint that each operator must act on 
two anti-parallel spins). Successive configurations main- 
taining the constraints are generated with three types of 
updates. 

In the first update — the "diagonal update" — the com- 
bined string P^L = (Pg) T Pa (where the transpose T of 
an operator sequence just corresponds to writing it in the 
reverse order, corresponding to acting with it on a bra 
state instead of a ket) containing 2m operators is tra- 
versed and each diagonal operator in it is updated (moved 
to a randomly selected bond), under the condition that 
it acts on anti-parallel spins. This step corresponds to 
changing the vertex break-up in the original world-line 
loop scheme*^ As in the SSE method^^ the constraints 
are checked by keeping the single state Z{p— 1), which is 
needed for moving a diagonal operator at location p in the 
string. This state is obtained by acting on the originally 
stored ket spin configuration Z(0) = Z r with the first p 
operators in the sequence. It is changed (by flipping two 
spins) whenever an off-diagonal operator is encountered 
in the course of traversing the positions p — 1, . . . , 2m. 
At the end of this procedure the stored bra state is ob- 
tained, Z(2m) = Zi : for a valid configuration. 

In a second updating stage — the loop update — a linked 
list of operator- vertices is first constructed. A vertex con- 
sists of the spin states "entering" and "exiting" an oper- 
ator, as shown in Fig. [U They connect, forming loops. 
The only difference with respect to the operator-loops in 
the SSE method is that a loop can now be connected to 
the ket or bra VB state, and the valence bonds consti- 
tute parts of such loops (replacing the periodic boundary 
conditions used at T > 0). To keep nonzero (indeed, con- 
stant) matrix elements of the operators H a b, all spins on a 
loop have to be flipped together, in the process changing 
also H a b(l) -H> H a b(2). Each loop is flipped with proba- 
bility 1/2. In practice, all loops are constructed, and the 
random decision of whether or not to flip a loop is made 
before the loop is constructed. Vertices in a loop that is 
not to be flipped are just flagged as visited, so that the 
same loop is not traversed more than once (i.e., a loop 
construction is always started from a vertex-leg that has 
not yet been visited). 

The reason for constructing all the clusters and flip- 
ping each with probability 1/2, instead of generating sin- 
gle clusters starting from random seed locations and flip- 
ping them with probability 1 (as in the classical Wolff 
method 2 -), is that the de facto loop structure is only 
changed when performing the diagonal updates. One 
would therefore potentially generate the same cluster sev- 
eral times, which would lead to lower efficiency compared 
to uniquely identifying all clusters and flipping each at 



most once. In principle one could modify the algorithm 
with combined diagonal and cluster updates, but this is 
more complicated and would probably not lead to im- 
provements in efficiency in most cases. 

A flipped loop including one or several VBs will cause 
spin flips in the stored spin configurations Z l or Z T . In 
the loop updating procedure we do not have to explic- 
itly keep track of any other spins than those in Z l and 
Z r . The four spins at the operators (the vertex legs) 
are irrelevant at the loop updating stage, because all the 
vertices automatically involve only operations on anti- 
parallel spins, both before and after a loop flip. For 
each vertex encountered when constructing a loop, we 
therefore simply have to change the operator-type index, 
1 ++ 2, in the list of operators (i.e., the same list P^L 
used in the diagonal update and to construct the linked 
vertex list). 

The third type of update — the state update — is iden- 
tical to the VB reconfigurations described in Sec. IIIII for 
the variational calculation. Normally one would use an 
amplitude-product state with coefficients (J9)), which en- 
ter in the weight (|2"o| . Reconfigurations of the bonds 
can be carried out with either two-bond or bond-loops 
moves, as explained in Sec. IIIII They only change the 
loop connections at the VB "end caps" . 



D. Measuring observables 

When measuring operator expectation values one can 
go back to a pure VB (=loop) representation, using the 
estimator (f2"3")l . This corresponds to summing over all 
loop orientations. Most quantities of interest can be ex- 
pressed in terms of the loops in the transposition graph 
corresponding to (V^(/3)|V r r (a))iSiS2i22i22i Note that these 
transposition-graph loops can also be obtained from the 
"space-time" loops constructed in the updates, by con- 
necting the sites (in practice just assigning a label, the 
loop number Aj) crossed by the same loop at the prop- 
agation midpoint (indicated by a dashed line in Fig. [4}. 
The " space-time" loops can also provide access to imagi- 
nary time correlation functions 2 - in the ground state (see 
section fVI A[) . Since there are no differences in the mea- 
surement procedures for equal time observables with re- 
spect to the original VB projector algorithm, we refer to 
the literature for this aspect of the simulations! 20 ! 21 ' 33 

In some applications, instead of measuring a ground 
state expectation value (0|A|0) one is interested in ma- 
trix elements of the form (H|A|0), where \R) is a refer- 
ence state, normally the Neel state in the z-component 
basis. This corresponds to sampling the wave function it- 
self (generating the basis states with probability propor- 
tional to the positive-definite wave-function coefficients) . 
The energy (including excitation energies in different mo- 
mentum sectors) can be computed like this , 20 ' 21 and also 
calculations of entanglement entropy can be formulated 
in this wayi 38 ' 39 ' 41 A mixed matrix element can also easily 
be sampled in the spin-bond basis. In this case the loops 



terminating on the state \R) should never be flipped, be- 
cause \R) is a single spin configurations (in the case of 
the Neel state — other reference states are also possible 
and would require other rules for the boundary loops). 



V. RESULTS 

As a demonstration of the efficiency of the methods, 
we present results for the sublattice magnetization M a of 
the 2D Heisenberg model. This quantity has been cal- 
culated in numerous previous studies, but the currently 
best published estimate, M a = 0.3070(3), obtained on 
the basis of T ~ QMC results for L up to 16, is already 
more than ten years oldi^ 2 - Recently, the density matrix 
renormalization group method was used to calculate M s 
on rectangular lattices with N « 200 sites, giving a re- 
sult consistent with the above value and with a similar 
precision^ 3 Results have also been obtained using finitc- 
T data and scaling forms that in principle allow simulta- 
neous T — > 0,L — > oo extrapolations. With L up to 160 
and 1/T up to 12, Ref. 44 reported M a = 0.30793(3). 
This is higher than (and well outside the error bars of) 
the T = results cited above. In order to resolve the dis- 
crepancy, it would be useful to have ground state results 
based on larger lattices. Here we consider L up to 256. 

Below we first discuss convergence aspects of the VB 
method, including the behavior with different trial states, 
and then present results and finite-size extrapolation of 
the sublattice magnetization. 



A. Variational calculations 

We first discuss the amplitude-product states used as 
trial states for the ground-state projection. The quality 
of the variationally optimized states [i.e., all amplitude 
h(x, y) were determined by variational Monte Carlo sim- 
ulations, as explained in Sec. IIIIj is illustrated in Fig. [5] 
for system sizes L up to 80. Results for up to L = 32 were 
previously presented in Ref. [24j — here we improve slightly 
on those results, thanks to the more efficient sampling 
procedures allowing for better statistics for the computed 
derivatives. The results are compared to converged re- 
sults of the QMC projector method (which can be con- 
sidered as exact to within small statistical errors that are 
not visible in the graphs) . The relative error of the varia- 
tional energy is < 0.1% for large systems. The sublattice 
magnetization falls on a a smooth curve in good agree- 
ment (better than 1%) with the projected data for L up 
to « 24. For larger systems the behavior becomes erratic, 
however, being higher or lower (outside the error bars) 
than the projected data in a seemingly random way. This 
can be explained as due to the energy becoming less sen- 
sitive to the long-range spin correlations for increasing L, 
i.e., there are states with significantly different sublattice 
magnetizations but energy expectation values that are 
the same to within the precision of the simulations. To 
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FIG. 5: (Color online) The energy (lower panel) and the 
squared sublattice magnetization (upper panel) of the vari- 
ational and ground-state projected states. 

obtain the correct best sublattice magnetization for large 
L (corresponding to the minimum energy determined to 
extreme precision) with the variational approach there- 
fore requires unreasonably long simulations (which is true 
in general in variational calculations; not just with the 
amplitude-product states used here). 



B. Convergence of the ground-state projection 

Turning now to results of the projector method, it is 
useful to test the convergence as a function of the projec- 
tion power to for different trial wave functions. Clearly, 
the preferred option is to use the best variational state 
available, but optimizing an amplitude-product state also 
takes some time (depending on how close to the energy 
minimum one strives), and, as we have seen above, for 
large systems it may not even be possible to find the 
truly optimal amplitudes. Fig. [5] shows the energy and 
the sublattice magnetization for L = 32 versus to/TV, ob- 
tained using trial states with amplitudes h(r) — l/r p , 
p = 2, 3, 4, without any optimization, as well as with am- 
plitudes obtained in two independent optimization runs. 
It is know n 24 ' 42 that the optimal amplitudes decay as 
1/r 3 asymptotically, but the short-bond amplitudes show 
deviations from this form. Indeed, the best convergence 
is seen for p = 3, but with optimized amplitudes the 
convergence is still much faster. Although the two opti- 
mized variational states have very similar energies, there 
are still clear differences in the convergence of the sub- 
lattice magnetization, related to the insensitivity of the 
variational energy to the long-distance spin correlations. 

In some cases the convergence of the sublattice magne- 
tization is non-monotonic (while the energy always has to 
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FIG. 6: (Color online) Convergence of the energy (lower 
panel) and the squared sublattice magnetization (upper 
panel) for L — 32 states states projected using different trial 
states; amplitude-product states with amplitudes h(r) = l/r p 
(p — 2, 3, 4) as well as with h(x, y) determined by minimizing 
the energy (in two independent optimizations, giving slightly 
different amplitudes). 



0.1050- 



0.1045- 



0.1040- 



0.1035 



0.1030, 



- ' ' ' 


i i i 


i ' i ' ' . 


- L = 64 


u 


i 


- 


- 


0.1250 


L = 20 


. 


- 


■ 




i 




: 


- 


0.1248 




. 


- 


* 




; \V**~ 


=« = •=*= A* 


■ 


-* 
_» 


0.1246 


1* 


- 


- 


* 


2 4 6 


8 10 


W 


^» »»--*-»--« 


.: 




■ 


) 5 


10 


15 20 


25 30 





m/N 



FIG. 7: (Color online) Convergence of the squared sublattice 
magnetization for L — 64 (L — 20 in the inset), using an 
optimized trial state. The dashed lines show the result ± 
error bar of SSE calculations (using loop updates) at very 
low temperatures (/3 = 8192 in the case of L = 64). 

converge monotonically), as illustrated in Fig. [3 The be- 
havior depends on details of the variationally optimized 
amplitudes; likely non-monotonicity can be traced to in- 
complete optimization. 



C. Extrapolation of the sublattice magnetization 

We now discuss large-scale calculations for the 2D 
Heisenberg model. We have calculated M 2 as well as the 
spin correlation function C(L/2,L/2), which equals M 2 



L 


Ml 


C(L/2,L/2) 


8 


0.177843(1) 


0.137595(2) 


10 


0.159372(2) 


0.128552(2) 


12 


0.147448(2) 


0.122586(2) 


14 


0.139153(2) 


0.118380(2) 


16 


0.133067(2) 


0.115263(2) 


18 


0.128412(2) 


0.112857(2) 


20 


0.124748(2) 


0.110954(2) 


24 


0.119350(2) 


0.108125(2) 


28 


0.115573(2) 


0.106126(2) 


32 


0.112782(2) 


0.104636(2) 


40 


0.108943(3) 


0.102571(3) 


48 


0.106431(3) 


0.101208(3) 


56 


0.104661(3) 


0.100239(3) 


64 


0.103345(3) 


0.099514(4) 


80 


0.101523(4) 


0.098501(4) 


96 


0.100325(5) 


0.097831(5) 


128 


0.098843(16) 


0.096990(17) 


192 


0.097371(11) 


0.096161(11) 


256 


0.096669(17) 


0.095765(16) 



TABLE I: Projector QMC results for the squared sublattice 
magnetization and the spin correlation function at maximal 
separation for several L x L lattices. The numbers within 
brackets indicate the statistical error (one standard deviation 
of the average) in the last digit. 

when L — ¥ oo, for lattices with L up to 256, making sure 
that the results are well converged to the ground state in 
all cases. The raw data are listed in Table, fl] The results 
are graphed versus 1/ L in Fig. [8j along with polynomial 
fitsii used to extrapolate to L — oo. The extrapolated 
M 2 and C(L/2, L/2) agree statistically and are stable 
with respect to the range of L included and the order of 
the polynomials. The statistics is slightly better for C 
and the polynomial needed to fit it is one order smaller 
than for Af s 2 . Based on C, we estimate M s = 0.30743(1), 
somewhat above the previous T = rcsultSi 12 ' 43 The er- 
ror bar is more than an order of magnitude smaller. The 
higher value from finite- T simulations— can be ruled out 
(differing by more than 15 of its error bars from our re- 
sult). This illustrates difficulties with unknown correc- 
tions to the (T, L) scaling forms. Extrapolating T = 
properties directly as a function of a single parameter 
(1/L) can in general be expected to be more reliable. In- 
deed, since the appearance of the (unpublished) original 
short version of the present article,— the results of Ref.|4J 
have been recalculated and revised 3 ^ and now agree with 
our results. 



VI. SUMMARY AND DISCUSSION 

We have shown how re-introducing the spins into the 
VB basis allows for much faster sampling of amplitude- 
product states and, in particular, very efficient ground- 




FIG. 8: (Color online) Finite-size scaling of the sublattice 
magnetization. The curves are polynomials fitted to 16 < 
L < 256 data (cubic for C and 4th-order for Ml ). The inset 
shows the deviation of the simulation results for C(L/2, L/2) 
from the corresponding fit. 

state projector simulations. For variational calculations 
it saves a factor up to system size N, and for projector 
QMC a factor up to the projection power m, bringing the 
total computational effort from 0(m 2 ) to 0(m) (where 
normally m 3> N, i.e., the improvements can be orders 
of magnitude). One striking and appealing aspect of the 
projector algorithm is its close similarity with state-of- 
the-art T > QMC loop methods, in particular the 
SSE variant^ but also the world-line approach (for which 
loop updates similar to those used here were originally 
developed^ 2 -). Essentially, the T — projector approach 
corresponds to "cutting open" the periodic imaginary- 
time boundary condition used at T > and terminating 
the resulting free loop ends with valence bonds (which act 
as continuations of the loops). It is therefore very easy to 
rewrite, e.g., an SSE program for T = projection. A fa- 
vorable aspect of the T = projector approach is that VB 
amplitude-product states often are very good trial states 
(as noted already a long time ago 23 i 25 ), which helps sub- 
stantially to achieve a fast convergence as a function of 
the power m of the projection operator H m . 



A. Discussion 

In view of the similarity with T > methods, it is 
worth asking how much faster a converged T = simu- 
lation is than a T > simulation carried out at T suf- 
ficiently low for obtaining ground state properties. To 
answer this question, we first need to compare directly 
the size of the configuration space of the two formula- 
tions. Consider the SSE method, which is based on sam- 
pling the Taylor expansion of e^^ H ^- The average of 
the expansion power n is given by (n) = /3\E\, where 
E is the total energy (ex N), which includes for each 
bond (in the case of the S = 1/2 Heisenberg model) 
an added constant —1/4 [exactly as in the singlet pro- 
jector operator in Eq. (|T5|) ]. Let us consider the 2D 
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FIG. 9: (Color online) Convergence of VB projector results 
for the energy (bottom panel) and the sublattice magnetiza- 
tion (top panel) for an L — 64 system versus the power m of 
the Hamiltonian used in the projection. The results at m/N 
are compared with T > results graphed versus the nor- 
malized inverse temperature 0.6/3 (corresponding to a similar 
computational effort as the projection scheme at m/N). 

Heisenberg model, where for large systems the energy 
per site is ss —0.67, which, including the constant —1/4, 
corresponds to \E\/N « 1.2 Thus, the average expan- 
sion power (n) sw 1.2N/3. In the projector approach, the 
number of operators in the sequence is 2m, and thus the 
computational efforts of the two methods should be com- 
parable if m/N ss 0.6/3 (since there is only a very small 
overhead in sampling in the projector scheme, related to 
the boundary VB states) . In Fig. [9] we show results of 
the two methods versus these normalized parameters for 
L = 64. A much faster convergence of the projector ap- 
proach can be seen, especially for the energy. It is clear 
that here the fact that the trial state is well optimized 
plays a big role — with a very poorly optimized trial state 
the T > approach may even converge faster. 

It should be noted that the loop estimators used in 
QMC calculations at T > are exactly equivalent to VB 
estimators^ such as Eq. Q, if T is low enough for a sim- 
ulation to sample only the singlet subspace, whereas at 
higher T some loops — those that change the total mag- 
netization when flipped — correspond to higher spin con- 
tributions. 

An important difference between the T — and T > 
approaches should be pointed out in this context: To im- 
prove the statistics in T > simulations one can take 
advantage of the periodic "time" boundaries. Averages 
of observables can be computed over all the propagated 
states. Whether or not this is a significant advantage in 
practice depends on the observable (its auto-correlation 
function as a function of imaginary time, which deter- 



mines the effective number of independent measurements 
generated). In some cases, the averaging can be carried 
out almost without overhead. At first sight it might ap- 
pear that no averaging of this kind can be done with the 
projector method, since the matrix elements of observ- 
ables are defined at the mid-point of the configuration 
(as indicated in Fig. Q|. One can of course shift the mea- 
surement point, and average over several such points, al- 
though then either the ket or the bra might not be equally 
well converged to the ground state. In practice, this kind 
of averaging may, however, still pay off, although we have 
not investigated it quantitatively. Note also that quanti- 
ties given in terms of Kubo integrals (susceptibilities) ex- 
plicitly require access to imaginary-time correlation func- 
tions, which can in principle be measured in the projector 
scheme by using the "space-time" loops constructed in 
the updates.— It would then be better to formulate using 
the Taylor expansion of e~P H instead of a fixed power 
(—H) m , and one would require longer projections than 
equal-time correlations, long enough for the imaginary- 
time correlator to have decayed below some small num- 
ber beyond which further contributions can be neglected. 
These issues are common to projector methods and have 
been investigated and discussed in detail in the context 
of fermionic auxiliary-field QMC calculationsi 36 i 37 



B. Related calculations and outlook 

The possibility to use the valence-bond basis in QMC 
simulations was first suggested more than 20 years agoi^ 3 . 
It is also implicitly the basis for the loop algorithm, 2 
Other aspects of simulations explicitly formulated in the 
VB basis became more widely used only recently, after 
its generic utility, including simplified access to many 
physical quantities, was pointed out . 20 ' 21 ' 33 We briefly 
list some works where the unique aspects of the VB ba- 
sis where exploited in such simulations, and where the 
improvements presented here can be expected to lead to 
further progress. 

Alet et al£& and, independently, Chhajlany et al£2- de- 
fined a VB entanglement entropy, which can be evaluated 
using ground-state projection (unlike other entropy def- 
initions, which normally cannot be directly evaluated). 
A slightly different variant of this entanglement measure 
was proposed by Lin and Sandvik, who also introduced 
another measure, the loop entanglement entropy, based 
on the loops in the transposition graph. 40 Hastings et al. 
recently showed that one of the standard definitions of 
entanglement entropy, the Renyi 1S2 entropy, in fact also 
can be evaluated with VB simulations. 41 These develop- 
ments will allow, e.g., tests of the "area law" of entan- 
glement entropy and deviations from it in a wide range 
of quantum spin systems. 

Beach et al. have extended the projector scheme for 
standard SU(2) spins to an often used representation of 
SU(iV) invariant models, including also the possibility 
of continuously varying N^- This extension is of interest 
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as many analytical theories are formulated as large- N ex- 
pansions, and it is useful to make direct contact between 
this approach and unbiased numerical calculations. Us- 
ing a similar approach, Tran et al. have carried out VB 
simulations of a chain of non-abelian SU(2)fc particles 
(where k is related to the central charge),— generalizing 
the standard 5=1/2 Heisenberg chain. One can also 
consider effective hamiltonians explicitly constructed in 
the VB basis42 

One of the first applications of the pure VB projector 
method was to the J-Q model; a potential realization of 
"deconfined" quantum criticality^ Since the appearance 
of the unpublished original short version of the present 
article,— the improved VB loop method has already been 
used to study SU(7V) versions of the J-Q model with 
N = 3 and 44£ Here it can be noted that the varia- 
tional amplitude-product states can also be extended to 
allow for valence-bond-solid order (by introducing factors 
favoring or disfavoring various local short-bond configu- 
rations), which is useful for studying J-Q models in the 
VBS phased 

The spin texture around a non-magnetic impurity in 
the J-Q model has been studied with the spin-bond sam- 
pling algorithm with an extra unpaired spin, 50 Triplet 
excitations of random antiferromagnetic clusters were in- 
vestigated with VB simulations including a triplet VB in 
the background of singlets by Wang and Sandvik^i a 
measurement which also can be carried out with two un- 
paired parallel spins. Measurements of triplet excitations 
for fixed momentum were discussed in Ref. l2ll 

Finally, we note that the loop-based projector scheme 



that we have presented here for isotropic spin systems 
can in principle be extended to anisotropic models as 
wel&£ (and also to more complex bosonic systems, as 
well as fermions in one dimensions). The difference will 
be that the variational trial states appropriate in most 
other cases would not have natural expressions in the sin- 
glet valence bond basis. As an example, consider the case 
of a variational wave function written just in terms of the 
z-components of the spins (e.g., a Jastrow state). There 
is then some weight (wave function coefficient) W({Sf}) 
for the trial state, and normally this weight will change 
when a loop terminating at the trial state is flipped. 
The loops in this case would be of the "directed" type^ 
and one could build in the weight-change coming from 
the trial state into the directed loop equations (which 
dictate the probabilities for the loop-building along dif- 
ferent paths through the vertices). This approach may 
be more efficient than the Green's function or diffusion 
Monte Carlo methods 52 normally used for ground-state 
projection of a variational trial state. 
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